%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib.pyplot as plt
import IPython.display as ipd
If we have an "analysis signal" $x[n]$ with $N$ samples, we define the Discrete Fourier Transform (DFT) by taking dot products with "probe signals," which are sines and cosines that go through an integer number of periods over the interval of $N$ samples.
def get_dft(x):
N = len(x)
t = np.linspace(0, 1, N+1)[0:N]
n_freqs = int(np.ceil(N/2))
cos_sums = np.zeros(n_freqs)
sin_sums = np.zeros(n_freqs)
for i in range(n_freqs):
c = np.cos(2*np.pi*i*t)
s = np.sin(2*np.pi*i*t)
cos_sums[i] = np.sum(c*x)/(N/2)
sin_sums[i] = np.sum(s*x)/(N/2)
return cos_sums, sin_sums
N = 100
n = np.arange(N)
x = 2*np.cos(2*np.pi*(3)*n/N)
x += np.sin(2*np.pi*(5)*(n-2)/N)
c, s = get_dft(x)
plt.figure(figsize=(12, 3))
plt.plot(x)
plt.title("Original Signal")
plt.figure(figsize=(8, 4))
plt.subplot(121)
plt.stem(c[0:10])
plt.title("Cosine Components")
plt.subplot(122)
plt.stem(s[0:10])
plt.title("Sine Components")
Text(0.5, 1.0, 'Sine Components')
The inverse Discrete Fourier Transform is going backwards; that is, given the result of doing dot products with all cosine and sine probes, go back and reconstruct the analysis signal. This can be done with the following formula and code
def get_idft(c, s):
"""
Given the cosine and sine components, reconstruct the signal x
in the time domain
"""
N = 2*len(c)
n = np.arange(N)
x = np.zeros(N)
n_freqs = len(c)
for k in range(n_freqs):
x += c[k]*np.cos(2*np.pi*(k/N)*n)
x += s[k]*np.sin(2*np.pi*(k/N)*n)
return x
We then check to see if our answer agrees with the signal we started with, and indeed it does!
x_recon = get_idft(c, s)
plt.figure(figsize=(12, 4))
plt.plot(x_recon)
plt.plot(x, linestyle='--')
plt.legend(["Reconstructed Signal", "Original Signal"])
<matplotlib.legend.Legend at 0x7f9af6111130>
from risset import *
def do_risset_faster(filename, tune_length, freqs_per_note, sr):
"""
Implement the faster version of Risset beats that aggregates
duplicate frequencies into a sine and cosine term, and then
uses the fast inverse Discrete Fourier Transform to reconstruct
the time signal
Parameters
----------
filename: string
Path to file with the tune
tune_length: float
Length, in seconds, of the tune
freqs_per_note: int
Number of frequencies to use for each note
sr: int
The sample rate of the entire piece
"""
ps, times = load_tune(filename, tune_length)
N = int(sr*tune_length)
c = np.zeros(N) # Analagous to dictionaries
s = np.zeros(N) # but using units of cycles per interval (instead of hz, cycles/sec)
for p, shift in zip(ps, times):
fhz = 440*(2**(p/12)) # In cycles/second (hz)
k_center = int(fhz*N/sr)
# Goal with risset beats to create a beat frequency that happens once
# over our tune interval
# Frequency k_center + 1 cycles/interval, difference is 1 cycle/interval
for dk in range(-freqs_per_note//2, freqs_per_note//2+1):
k = k_center + dk
f = k*sr/N
c[k] += np.cos(2*np.pi*f*shift)
s[k] += np.sin(2*np.pi*f*shift)
## Setting up format to use the "Fast Fourier Transform," which is an O(N log N) algorithm
## whereas the "brute force" version we implemented is O(N^2)
M = N//2
if N%2 == 1:
M = (N-1)//2
else:
M = (N-2)//2
c[-M::] = c[1:M+1][::-1]
s[-M::] = -s[1:M+1][::-1]
return np.real(np.fft.ifft(c - 1j*s))
import time
tune_length = 70
freqs_per_note = 500
sr = 44100
tic = time.time()
x = do_risset_faster("Tunes/wanna.txt", tune_length, freqs_per_note, sr)
print("Elapsed Time: {:.3f}".format(time.time()-tic))
ipd.Audio(x, rate=sr)
Elapsed Time: 0.241
Let's look at the difference
N = 44100*70
# N^2 / N log2(N), N/log2(N)
N/np.log2(N)
143196.6024188806